5 ' L 


JPO, 1999 


Ocean Turbulence I: one-point closure model 
Momentum and heat vertical diffusivities 


V.M. Canuto 1 ’ 2 , A. Howard 1 , Y. Cheng 1 and M. S. Dubovikov 1 
! NASA Goddard Institute for Space Studies 
2880 Broadway, N.Y., N.Y., 10025 
2 Dept. Applied Physics, Columbia University, NY, NY, 10027 


1 



Abstract 


Since the early forties, one-point turbulence closure models have been the canonical 

tools used to describe turbulent flows in many fields. In geophysics, Mellor and Yamada 

applied such models using the 1980 state— of-the art. Since then, no improvements were 

introduced to alleviate two major difficulties: 1) closure of the pressure correlations, which 

affects the correct determination of the critical Richardson number Ri above which 

cr 

turbulent mixing is no longer possible and 2) the need to express the non-local third-order 
moments (TOM) in terms of lower order moments rather than via the down— gradient 
approximation as done thus far, since the latter seriously underestimates the TOMs. Since 
1) and 2) are still being dealt with adjustable parameters which weaken the credibility of 
the models, alternative models, not based on turbulence modeling, have been suggested. 

The aim of this paper is to show that new information, partly derived from the newest 
2-point closure model discussed in paper III, can be used to solve these shortcomings. The 
new one-point closure model, which in its simplest form is algebraic and thus simple to 
implement, is first shown to reproduce a variety of data. Then, it is used in a O-GCM 
where it reproduces well a large variety of ocean data. 

While phenomenological models are specifically tuned to ocean turbulence, the present 
model is not. It is first tested against laboratory data on stably stratified flows and then 
used in an O-GCM. It is more general, more predictive and more resilient, e.g., it can 
incorporate phenomena like wave-breaking at the surface, salinity diffusivity, non-locality, 
etc. One important feature that naturally comes out of the new model is that the predicted 
Richardson critical value Ri £r is Ri ~1 in agreement with both LES and empirical 
evidence while all previous models predicted Ri cr ~0.2 which led to a considerable 
underestimate of the extent of turbulent mixing and thus to an incorrect mixed layer 
depth. The predicted temperature and salinity profiles (vs. depth) are presented and 
compared with those of the KPP model and Levitus data. 
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I. Introduction 


One-point closure turbulence models have a long tradition that begun with the 
pioneering work of Chou (1940, 1945) that was subsequently extended and improved by 
Lumley and Khajeh-Nouri (1974), Launder, Reece and Rodi (1975), Pope (1975), Lumley 
(1978), Rodi (1984) and Shih and Shabbir (1992). For a review, see Gatski et al., (1991). In 
a pioneering work, Mellor and Yamada (1982, MY; Mellor, 1989) applied the 1980 
state-of-the-art 1-point closure models to geophysical problems. The original MY model 
is known for its successes and its failures: 1) it does not parameterize the pressure 
correlations correctly. This has two consequences: it predicts a critical Richardson number 
Ri cr =0.2 (above which turbulence can no longer be maintained) that is much smaller than 
Ri cr ~l required by empirical evidence and LES (large eddy simulations, Wang et al, 1996), 
which means that the YM model underestimates the extent of turbulent mixing and given 
in correct mixed layer depth, second, the model does not fit the PBL surface data, 2) it 
treats the third-order moments (TOM) with the down-gradient approximation (DGA) 
which is now known to severely underpredict the true TOM. And yet, even the most recent 
applications of the MY model to ocean turbulence (e.g., Rosati and Miyakoda, 1989; 
Galperin et al., 1989; Gaspar et al., 1990; Blanke and Delecluse, 1992; Baum and Caponi, 
1992; Ma et al., 1994; Kantha and Clayson, 1994; Burchard and Buamert, 1995; D'Alessio 
et al., 1998), still use the original model with changes that do not solve the difficulties 1) 
and 2). It is fair to say that, unless and until such shortcomings are solved, the credibility 
of the model is at stake. We have solved 1) and 2) as follows. Pressure correlations 
originate from the Navier-Stokes equations for the turbulent velocity field which contain 
the pressure gradients. For example, in deriving the equations for the Reynolds stresses and 
the heat fluxes, vectors and tensors appear which represent the correlation of fluctuating 
temperature and velocities with p r. 

n?E^ iy^ + 3^ (i) 


3 



where p. and 0 are the fluctuating components of the pressure, velocity and temperature 
fields. Since the fluctuating pressure p does not satisfy the hydrostatic equilibrium 
equation, these higher-order terms must be modeled in terms of second-order terms. As 
discussed in Appendix A, the most complete expressions are (for sake of simplicity, we 
leave out the numerical coefficients; all the tensors are defined in Appendix A) 

( 2 ) 


V 2 >Vf K V B ij-V z ij 

n f " T P «¥ +A i ?! -( s ij + 5 v ij)¥ 


( 3 ) 


The first terms in both 2) and 3) represent the return to isotropy or Rotta terms; the 
second and third terms in (2) are the contributions due to shear (S—) and buoyancy (B- j); 
the fourth term is due to the non-isotropic contribution due to shear while the last term is 
entirely due to the vorticity \B. None of the second-order closure models employed thus 
far in 0— GCM has employed the full form (2-3) and that is why these models often 
underestimate the true value of R' cr - For example, even the rather complete analysis by 
Kantha and Clayson (1994) adopts the simplified expressions 


Ilf = r' 1 / ) uT5+ A.# 2 — S : .u7# 


n.. = 2r~ ! b.. - -KS-. 

ij pv ij 5 ij 


p0 u i 


»J J 


(4a) 


which neglects the last three terms in IB. The model of D'Alessio et al. (1998) adopts an 
even simpler expression: 


n ? =r pi>¥ +A i 5! 


(4b) 


which also neglects terms in IB which we have found instrumental in the determination of 
Ri Cf . We have adopted the full forms (2)— (3). 


The second problem concerns the pressure-temperature and the pressure-velocity 


time scales r t . Specifically, the question arises as how to evaluate the ratios 

V /t ' v /r (4c) 

assuming that the dynamical time scale r=2Kf _1 is known from solving the dynamic 
equations for K and t. All one-point closure formalisms thus far have been unable to 
determine (4c) and empirical constants were thus introduced to represent r «/r and y /r. 
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These adjustable parameters weaken the credibility of the model. The two— point closure 
model discussed and employed in paper I has allowed us to determine such time scales 


whose values are no longer free. They are given by: 

J3v 2 T p0 1/, . _jx_i t 9 , AAX 

T =5> T =5< 1 + V> - T =a t < 4d ) 

where a ^ = 0.72 and where we have added the temperature variance (e# 2 ) dissipation time 
scale t q. An immediate bonus of these improvements is that the new model matches the 
measured surface PBL data rather well: 


Table I 



u 2 /2K 

v 2 /2K 

w 2 /2K 

Ri(cr) 

PBL data 

0.56 

0.27 

0.18 

0.5 

MY Model 

0.56 

0.22 

0.22 

0.2 

New model 

0.56 

0.27 

0.18 

0.5 


Next, we discuss the problem of the third-order moments TOM which represent 
non— locality. The prototype is the flux of turbulent kinetic energy: 

F(K) = iq 2 w, K =iq 2 (5a) 

All ocean turbulence models used thus far either disregard the TOM altogether (local 
model), or use the DGA which assumes that the kinetic energy flux is along the gradient of 
K: 

F ( K )=- K m* < 5b > 

However, as shown by Moeng and Wyngaard (1989), the DGA severely underestimates the 

true value of the TOM and thus, it should be avoided. All ocean turbulence models that 

employ a prognostic equation for K (Rosati and Miyakoda, 1989; Gaspar et al., 1990; 

Blanke and Delecluse, Ma et al., 1994) have adopted the DGA, (5b). In addition, the 

diffusivity K m is taken proportional to w t however, in the case of stratification, K m is a 

combination of both the momentum and the heat fluxes, 
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(5c) 


K m = aw/ + bw# 

and since for stably stratification, W9 is negative, the second term decreases the diffusi vity . 
Once the down-gradient approximation is abandoned, the only alternative is to consider 
the full dynamic equations for the third-order moments (Appendix B) and try to solve 
them. In the ID case, this was done and the results are as follows (Canuto et al. , 1994): 

q 2 w = E f- w 0 + E f- w 2 + E + E |V 
n 1 dz 2 dz i dz 4 


w 3 = B I-W0+B t-w 2 + B tf 2 + B Or- q 2 
i dz 2 dz 3 dz 4 az n 


w 2 6> = a|-w0+A |-w 2 + A ^-flP + A ^q 2 
\dz 2 dz 3 dz 4 dz M 


W0 2 = C ^-W , v, n 

1 dz 2 dz 3 dz 


2 + c i^ + C.^q 2 


=D. fj w« + D. L w 2 + D„ + D. |;q 2 


dz 


3 dz 


dz 


0 s = F |-wf?+F 4w 2 + F l-^ + F f-q 2 
1 dz 2 dz 3 dz 4 dz ^ 


(5d) 


Eqs.(5d) exhibit a remarkable symmetry: all the third-order moments have the same 
structure since they are linear combinations of the derivatives of all the second-order 
moments. Second, all the turbulent diffusi vities A-F exhibit the same general structure 
(5c). The DGA corresponds to taking 


A = 0 , B =0, C =0 

2i3?4 1}3 i4 1)2?4 

D = 0 , E = 0 , F =0 (5e) 

2,3,4 1,2,3 1,2,4 

In Canuto et al. (1994) the new TOM (5d) were tested against large eddy simulation data 
for the PBL. The results were quite good. Eqs.(5d) are the most complete representation of 
the TOM presently available. 

Finally, we discuss the dissipation length scale. The dissipation rate c of the turbulent 
kinetic energy K is one of the most crucial and most difficult variables to model. Although 
an exact differential equation for e was derived long ago, it has been used only very seldom 
(e.g., Burchard and Buamert, 1995). Most frequently, use is made of the local form: 

t = l A K 3 / 2 (5f) 
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The problem then reduces to the evaluation of the dissipation length scale l. Until recently, 
several empirical expressions for l have been proposed which often seem to contradict each 
other. Among these formulae, the most widely used are the ones due to Deardorff (1980) 
which includes only buoyancy, and by Hunt et al. (1988) which includes only shear. Neither 
expression is valid in the limit of small buoyancy or shear. Starting from a theoretical 
formulation of the buoyancy spectrum in stably stratified flow due to Lumley (1964) and 
Weinstock (1978), a new expression was derived for f that: a) includes both buoyancy and 
shear, b) encompasses all previous models of £, c) behaves well for small buoyancy and/or 
shear, and d) explains the non— monotonic behavior of l found by LES. 

In conclusion, we show that a one— point closure model free of the uncertainties of the 
past can be constructed. Tests of the model are presented before being used in an O-GCM, 
the results of which are also presented. On the basis of its overall performance, the model 
will be further expanded to include salinity (paper II). 


III. Turbulence model 

We begin by recalling that in an ocean GCM the equations for the mean velocity and 


mean temperature entail the Reynolds stresses UjUj and the heat fluxes In 


r if Vj’ 


h r¥ 


(6a) 


where Uj and 6 are the fluctuating components of the velocity and temperature fields. The 
well established procedure of deriving the dynamic equations for the variables (la) entails 
other turbulence variables, specifically (u and \ are the kinematic viscosity and thermal 
conductivity): 

K e jr.. (turbulent kinetic energy) (6b) 


V 2 = (temperature variance) 
i=v{ w f (dissipation rates of K) 


(6c) 

(6d) 
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~njj — 

t g = Xv^ ) (dissipation rates of F) (6e) 

1 

There are a total of 12 independent variables, each of which satisfies a separate dynamic 
equation, and each of which depends on the others so that it is not possible to derive the 
heat flux without solving the equation for the Reynolds stresses and viceversa. Since the 
derivation of the dynamic equations for the turbulence variables (6a— e) has been presented 
elsewhere (e.g., Gatski et al. , 1991; Canuto, 1994), we shall only cite the results. 


a) Turbulent kinetic energy K: 

PY + D f (K) = P s + P b - ( 

where the production terms P's due to shear and buoyancy are defined as 

P,e A.h., 
b 1 ■’ 


P S 5 - T ij S ij- 


i r 


(7a) 

(7b) 

(7c) 


while the divergence of flux of turbulent kinetic energy R(K) is given by 

DfWiJjiyn), Fj(K) = ^ 

Finally, AjEg.a^ with gj=(0,0,g) and a T =-p A dp/ d T. The shear is defined in Eq.(A.6) 

b) Temperature variance (potential energy) (P: 

£}2 _ 

(8a) 

P 0 = ~ h ifxj’ D f(i^)= ( 8b ) 

P^ and Dj- represent the rates of production and transport of the divergence of the flux of 
temperature variance and ty represents the rate of dissipation of i# 2 , Eq.(6e). 


m ^ 2 + D f^ ~ ? e~ c e + ** dxf 2 


o 

c) Reynolds stresses, b-. = r*. K&-: 

ij ij o 1J 


D 


DtVVV - -i5 K V 2r pl b ij + ^V 1 -“ 1 > s ij-< 1 -“ J ) z ij 


(9a) 


The traceless tensors B, £ and Z are defined as 


B..= A-h. + A.h- 4(5.. A, h, 
ij i J j i 3 ij k k 


(9b) 
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V s ikV s jk b ik-sWM 

z ij = v ik b kj +v jk b ik 
Shear S— and vorticity V- are defined as: 

s ij = *(££ + Up’ v ij = - l¥f 

J ^ J 


The diffusion term is defined as 


°f ( y - i ( u i u j u k~3 & ij^ u k) 


(9c) 

(9d) 

(9e) 

(9f) 


c) Heat flux: hj=m7J 


D, h i +I W - - r ijW. - tJ-7« 3 > h j S ij - ^tWu + <-'-V\ V - + ^ h i 

J ^ 

(10a) 

where 


<9T 


d 2 


d) Dissipation rate e q 

e) Dissipation rate c. 


D f (h i>= obclvi® + 


‘e = T e v 


gf + D f (r) = f K-Hc i P s + C 3 P b )-c/K- 


D f (f ) “ ijy 


(10b) 

( 11 ) 

(12a) 

(12b) 


Since we have assumed an algebraic for the number of dynamic equations is now 11. 


IV. Diffusion terms D^ 

The third-order moments are discussed in Appendix B. 


V. Algebraic Models 

The above model is complete and allows us to compute the turbulent variables of 
interest. In its full dynamical form, it is obviously rather complex since it entails 11 
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differential equations. In this section, we present a model in which we retain only two of 
the eleven differential equations, those for K and t. The procedure, known as ARSM. 
algebraic Reynolds stress models , neglects the D/Dt and the diffusion terms. After some 
algebra, we obtain: 


Reynolds Stresses: 


b.. = “TxKr S. ; + t B..— *r (l-a)E.-ir (1-a )Z.. (13) 

ij 15 pv ij 2 "5 pv ij 2 pv v R ij 2 pv v 2 B v ' 

Since the last two terms depend on b^ itself, Eq.(13)is a system of linear algebraic 
equations. 

Heat fluxes: 


A ik h k = -(K h)ii 


&I 


h'ij <7x ; 


J 


where the turbulent heat conductivity tensor is given by 


( K h)ii = ^ b ii + 


'B 


B 


B 


The dimensionless tensor A^ is defined as: 


A.. = - 
B t. 


d T 


P# 


*ij + ( 1_ V 7 Vi35T j + (1_ ? a 3 )rS ij +(1 ^ a 3 )rV ij 


(14a) 


(14b) 


(14c) 


The fact that the rhs of (14b) depends on b- j which in turn depends on In via B- j makes the 
analytic solution in the 3D case somewhat cumbersome though manageable with symbolic 
algebra. In order to homogenize the notation, we introduce the dimensionless constants: 

j * 

A=-^ v , A^= j 5 A, A 2 = i(l— a^A, X^= i(l-a 2 )A, A^= ±/? & A 

J 1 

A = l , A = 1- la , A = 1- la , A = ( 1-7 )— 

5 rtf 6 4 3 ’ 7 4 3 ’ 8 v \ T 

Thus, Eqs.(13) and (14c) can be w'ritten in the compact form: 


b.. = - A rKS.. + A rB..- A rE.. - A rZ.. 
B 1 B 4 B 2 B 3 B 

A-. = A &. + A ^X-M- + A rS i; + X rV-. 
B 5 B 8 idxj e B 7 B 


(14d) 

(14e) 

(14f) 


VI. Vertical Diffusivities 
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(15a) 


Here, we present the analytic solution of Eqs. (13) — (14) for the following case: 


ffi dn c 
~dx. - Tz ’ i 3 ’ 

The shear and vorticity acquire the form: 

r 0 0 d\J/dz 
S..=i 0 0 dV/dz , 

1J L dU/dz dW/dz 0 J 


q = [u(z), v(z), 0] 



0 0 dU/dz-, 

0 0 dV/dz 

t-d\J / dz - dV / dz 0 J 


(15b) 


Using symbolic algebra, one can solve Eqs.(14a,b,e, f). The results are as follows: 

Reynolds stresses: 


uw = - K 


dlj 

^mcfz' 

<9V 

'mcfz 

— / 1 . , j/ d\] dV 

uv- (A +A )rK m ^-^- 


vw 


-K 


Time scale: 



Turbulent kinetic energies: 


(16a) 

(16b) 

(17) 



^ = ! k + Wit v 3 v<!f> 2 - + kv N2 

(18a) 


* = J K + 5 TK ml<\+ 3 V<!f> 2 - 2A 2 ( ^ )2] + 3 A „ V" 2 

(18b) 


^=! K+ j(V 3 V K m rE2 -! A 4 v N2 

(18c) 

Mean shear: 

s2 * (|j) 2 +(|t) 2 

(18d) 

Heat fluxes: 

= ■ K hfz 

(19a) 


^ = - A ;'I K m + ' (A 6 + V K h! r H!F 

(19b) 


“ A 5 '[ K m + ’ (A 6 +A 7 )K hl T lz 1? 

(19c) 

Turbulent momentum and heat diffusivities: 



K m= K h = 2S h7 

(20a) 


As discussed in I, the dimensionless structure functions S m ^ are instrumental in 
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being able to show that the model satisfies a well known feature of stably stratified flows 
(Webster, 1964), namely that the turbulent Prandtl number 


K m 


is an increasing function of the Richardson number Ri: 

3Ki <r T (Ri) >0, Ri = ^f 2 ^ 
The function a (Ri) vs. Ri will be discussed in sec. XVIII. 


(20b) 


(20c) 


VII. Structure Functions 
Dimensionless functions S 


m,h' 


DS m = s + s (rN) 2 + s (rE) 
m o i 2 


Plots of S 


DS h = s +s (r N ) 2 + s (rE) 2 

11 4 5 6 

D = d + d (rN) 2 + d (rE) 2 + d (rN) 4 + d (r’NE) 2 + d (rE) 4 

0 1 2 3 4 5 

vs. Ri will be exhibited in sec. XVIII. 


m,h 

Dimensionless variables s^: 


|a A 2 

^ 1 5 


s e - A (A +A ) + 2A A (A - ,A -A ) + h A A 

1 4 6 7 4 5 1 2 3 ^ 1 5 £ 


S = - (A 2 ~A 2 ), S = 2A , 


-2X 


s = ?A (3A 2 -A 2 ) -fAA (3A -A ) + ?A (A -A ) 

6 1 5 3 1 5 1 3 2 7 4 i v 6 7 


Dimensionless variables d^: 


d = 3A 2 

0 5 


d = A (7A + 3A ), d e A 2 (3A 2 -A 2 ) - ?(A 2 -A 2 ) 

1 5 V 4 2 5' 3 1 4' 6 T 

d = A(4A +3A), d e}(A 2 - 3A 2 ) (A 2 -A 2 ) 

3 4 V 4 8 ; 5 4 v 2 3 y V 6 7* 

d = A[AA-3AA-A (A 2 -A 2 )] + A A (3A 2 -A 2 ) 

4 4 L 2 6 3 7 5 V 2 3 /J 5 8 3 2' 


(21a) 

(21b) 

(21c) 


(22a) 

(22b) 

(23a) 

(23b) 

(24a) 

(24b) 

(24c) 

(24d) 


As a way of comparison with recent work, we notice that the expressions for uw and w d 
employed by D'Alessio et al. (1998, see equations 6-12) are equivalent to taking E-*0, that 


is, their diffusivities do not depend on the source (shear). In that case our S m would 
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become a constant 


S m = 0.09 (24e) 

As we shall show, in our model this value corresponds to Ri<0, whereas one expects that 
S m decreases with Ri for the reasons given above. 

Since in the above relations K and e have remained unspecified, one can choose 
different models. 

VIII. First non-local model. 

In this case K and c are taken to satisfy the two non-local equations: 

DT + D(K) = K m £ 2 - K h Ng - < (25a) 

nf + D(e) = <K-‘(C K m S? - C 3 K h Ng) - c/K-' (25b) 

This type of model was adopted by Baum and Caponi (1992), Kantha and Calyson (1994) 
and Burchard and Baumert (1995). 

IX. Second non-local model. 

In this case, one retains only one prognostic equation, the one for the turbulent kinetic 
energy K, Eq.(25a), while Eq.(25b) is taken in its local limit 

f = A-'K 3 / 2 (25c) 

where the mixing length A has to be specified. Models of this type with down-gradient 
approximation for the diffusion term were considered by Rosati and Miyakoda (1989), 
Gaspar et al. (1990), Blanke and Delecluse (1993) and Ma et al. (1994). 

X. Third non-local model 

In this model, Eq.(25a) is taken in its local form by neglecting the left hand side 
altogether while one retains Eq.(26) so as to avoid the need to introduce a mixing length l. 
The local limit of (25a) is physically equivalent to assuming production=dissipation, 

< = K m E 2 - K h N 2 (26a) 
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or, using the relation r=2Kr 1 , 


(r£)’S m -(Ny) 2 S h = 2 

Using Eqs.(22a— c), Eq.(26b) becomes an equation for the dimensionless variable 

(rE) 2 M 

The result is: 


A tj? + Bip + C = 0 

A e (s +2d )Ri 2 - (s -s -2d )Ri - s + 2d 

5 3 M 6 r 2 5 

B = (s +2d )Ri - s + 2d , C = 2d 
4 r o 2 o 

The function ty vs. Ri will be discussed in sec. XVIII. 


(26b) 

(26c) 

(26d) 

(26e) 


XI. Fully local model 

In this case, both Eqs.(25) are taken to be local; K can then be expressed as 

£ = t(r\ K = 4A 2 E 2 (27) 

K o 0 

Once a model for A is provided, the model becomes fully algebraic. 

XII. The critical Richardson number 

The critical Richardson number Ri cf is defined as the value of Ri above which 
turbulent mixing is no longer possible due to the action of stable stratification. While linear 
analysis yields the result Ri cr =l/4 (Maslowe, 1981), it has been found (Martin, 1985) that 
a reliable prediction of the mixed layer depth can only be achieved if Ri ~1. Early 
laboratory data by Taylor (cited in Monin and Yaglom 1971) indicate that turbulent 
exchange exists even when Ri> 1 . In addition, recent LES studies (Wang et al., 1997) also 
show that turbulence exists up to Ri~l. This is therefore a critical test of any turbulence 
model. As discussed earlier, the original MY model gives an Ri cr =0.2 which is even smaller 
than the result of the linear stability analysis, Ri cr =l/4. We shall define Ri cr as the value 
at which the kinetic energy vanishes: 
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(28a) 


K(Ri cr )=0 

Since ^K' 1 , Eq.(26d) implies that in the limit K-*0, A(Ri cr )=0 which in turn gives 

Ri = (2c ) _1 [- c + (c 2 - 4c c )*] 
c v 1 2 2 1 3 J 

c = s + 2d , c = — s +s ■(■2d , c — — s ■+■ 2d 

15 32 1643 2 5 

Using the model constants discussed below, 

Ri £ = 0.846 (model A) 

Ri cr = 1.03 (model B) 


in agreement with the data referred to earlier. 


(28b) 


(28c) 


XIII. Realizability Conditions 

Since by definition (r£) 2 >0, from Eq.(26b) it follows that the minimum value of rN is 

( rN >mi„ = - 2S h' < 29a > 

or, using (21b) and the model constants, we derive: 

( rN >min = “ 12 - 27 (29b) 

2 

The negative value of ( r N) m j n means that it occurs in the unstable region. In the stable 
region, (rN) 2 is limited by stable stratification as follows (Deardorff, 1980) 

< TN )max = 466 < 29c > 

On the other hand, by arguing that an increase of shear should noi result in a decrease of 

-uw, one can derive an expression for the maximum value for (rE) 2 : 

< TE >max ~ l d 0 + d 1 ( rN ) 2 + d 3 ( rN ) 4 ][ <l 2 + d < ( , - N ) 2 ]' 1 
In applying the model, we shall use the following realizability conditions: 


(29d) 


(rN) 2 = (rN) 2 . 
v ’ v 2 mm 

\n 2 


if(rN) 2 <(rN)^ n 


( rN ) - ( rN W if ( 7 ’^) 2> ( rN ) max 


) 2 < 

) 2; 

\ 2 > 


2 

rri 

2 

m 

2 


( rE ) 2 =( rE )max if ( rE > 2 X rE )max 


(Me) 


XIV. The representation K^=7cN~ 2 

Several authors (Thorpe, 1977; Osborne, 1980; Oakey, 1982; Mourn, 1989; Davis et al., 
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1994; Toole et al., 1994; Wuest et al., 1996) have represented the turbulent heat diffusivity 
in the form: 

K h = 7^2 (30a) 

The derived values of 7 are: 


0.12<7<0.48 


Using (26a), we obtain: 


R. 


7 = 


f 

THr 


Rj. = Ri a^, 1 


(30b) 


(30c) 


where use must be made of (20b). The function 7(Ri) vs. Ri is discussed in sec. XVIII. The 
value 7=1/4 (e.g., Oakey 1982; Toole et al., 1994) occurs at Ri=0.25. In using (30a), care 
must be exercised here since in general the heat diffusivity is not the same as the "mass 


diffusivity" defined via the relation: 




It is easy to derive the relation 


K. 


R = — 
p a 


s dS/dz 


(31a) 


(31b) 


% = K h( 1 -<V (1 -V ’ 

where K s is the salt diffusivity and R^ is the Turner number (n g is the haline contraction 

r 

k 

P 


s ~ " "" * “ ~~p v “s 

coefficient). As in the present case, K g =K^, it follows that K^K • 


XV. Breaking Waves 

It has been known for a number of years (Kitaigorodskii et al., 1983; Gargett, 1989; 
Agrawal et al., 1992; Anis and Mourn, 1994; Drennan et al., 1996; Terray et al., 1996, 1997; 
Burgers, 1997; Skyllingstad et al., 1999) that in the ocean mixed layer beneath actively 
breaking waves the dissipation rate c of turbulent kinetic energy K is even two orders of 
magnitude larger than the one corresponding to the "law of the wall": 

f wall = U *( KZ )' 1 ’ u * = (^ w )* ( 32a ) 

which can be derived from (27). Here, where «~0.4 is the von Karman constant, z is the 

distance from the "wall", r is the surface wind stress and p w is the density of seawater. To 
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incorporate the wave breaking phenomenon into a turbulence model, one needs to solve the 
dynamic equation for the turbulent kinetic energy, Eq.(25). This in turn requires that a 
boundary condition to specify K at z=0. Craig and Banner (1994) first proposed the 
expression (K=^q 2 ) 

q | z _ 0 = (12c' 1 ) 1 / 6 u*/ 3 c 1 / 3 (32b) 

where c is the "effective phase speed" (Terray et al., 1996) that enters the rate of energy 
input to the waves from the wind, F=r c/fl ~u 2 CEmi* , where r =p u* =p u* 
(a:atmosphere, wrwater). Craig and Banner (1994; their equation 27) use a rather than c. 
The coefficient c^=0.18. In Appendix C we present a new and simpler derivation of (32b) 
and in sec. XVII I we discuss the results of our model with breaking waves. 

XVI. Modeling the length scale A 

We begin by writing A that appears in (25c) as 

A = c^Af(N,E) (33a) 

and require that f(0,0)=l. Here, A is the size of the largest eddy, k Q A=7r while f(N.E) 
represents the distortion of the Kolmogorov spectrum due to stratification and shear. 
Though a satisfactory theory of the energy spectrum E(k) under shear and stable 
stratification is still not available, there is some consensus (Gargett et al., 1981) about the 
overall shape of E(k): 

I ; E(k)=(eN)’k" 2 , Fr<l, Ri>l 

II: E(k)=cN 2 k- 3 , Fr=l, Ri~l 

HI : E(k)=Kof 2 / 3 k" 5 / 3 , Fr>l, Ri<l (33b) 

where Fr is the Froude number 

K * 

Fr = XN ( 33c ) 

For the ocean c»10 (Gargett et al., 1981; Gargett 1989), while atmospheric data suggest 

cslOO (Dewan, 1979) so that the interval between regimes I and II may be one or two 
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decades. Regions II and III coincide at a length scale given by the Ozmidov scale, 

A 0 = F = (r 3 )’ (33d ) 

Deardorff (1980), Hunt et al. (1988) and Fernando and Hunt (1996) have suggested the 
following phenomenological relations: 


K 


f(N,0) = 0.76 f(0,E) = 2.76 ^ (33e) 

Deardorff's result, the first of (33e), comes from substituting Eq.(27) in the Ozmidov scale. 
Neither expression satisfies the condition f(0,0)=l and thus they have a limited validity. 


Furthermore, f(N,0) implies that 
which in turn implies that 


~ AK'^ ~ N' 1 , rN ^constant 


(33f) 


0 = (rE) 2 ~ Ri' 1 (33g) 

while (28d) implies the opposite, namely a 0 that increases with Ri. Similar remarks are 
valid for the second of (33e). There have been several attempts to construct a model for 
E(k,N,E). Lumley's (1964) model implies a function f(N,E) of the form (Cheng and 
Canuto, 1994): 

f(N,E) = [1 - c Fr 2 (l - c ^Ri- 1 )]' 3 / 2 (34a) 

where cr = o- (R.i). The two constants are given by 

c = (27r 2 3’) 4 Ko 3 / 2 , C 2 =0.4-0.6 (34b) 

As discussed in Cheng and Canuto (1994), Lumley's formulation is valid only for small 
levels of stratification. Weinstock (1978) suggested an improvement of Lumley's model and 
the corresponding form of f(N,S) is (Cheng and Canuto 1994): 

f(N,E) = [(1-A 2 B 2 )’ - AC]- 3 (34c) 


where 


A = a (1 - c <7 Ri 4 ), B" 1 = a +Fr 2 , C 4 = a +Fr 2 

1 2 T n 2 3 


a 


= g.eslO^Ko 3 / 2 , a = |l0' 2 , a = i, 10- 2 


(34d) 


FTiU 5 Cl — 7) 

1 2 ^ 3 ^ 

Recently, Cheng and Canuto (1994) improved on both Lumley and Weinstock's model and 
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checked the result against LES data (see their Fig.5). The function f(N,£) is given by: 

f(N,£) = [1 + aFr-^l+bFr 4 / 3 )- 1 ]- 0 (34e) 

where 

R f <c : 
i 2 

a = 2(37r 2 )- 1 (n-l), b = 0.12 (0 - 1 — c= \ (34f) 

Rr> c : 

I 2 

a = 4(57r 2 ) 4 f2, b=0, c = ^ (1— O' 1 ) 

0= 1 +| ? Ko 3/2 ((7 T Ri- 1 -l) (34g) 

Canuto and Cheng (1997) have further improved on their original model and their new 
expression for f(N.E) is Equation (6a) of that paper. In the case in which one considers only 
shear, f(0,E) simplifies considerably: 

2f' 2 / 3 = 1 + pSh 2 + (1 + 2pSh 2 - ^p 2 Sh 4 )’ (34h) 

where 

Sh = AEK'% p=i(2^ 2 3’)- 1 Ko 3 / 2 (34i) 

In the simplest model, A is determined using the Deardorff-Blackadar formula: 

A = 2' 3 / 2 B L B^=24.7, 

l = min(i^, ^ ) (34j) 

i — Kit (1 4-kz)' 1 , l = 0.17H (34k) 

1 0 0 0 

where iq 2 =K is the turbulent kinetic energy, N is the Brunt— Vaisala frequency, k= 0.4 is 
the von Karman constant and H is the mixed layer depth. When used within the NCAR 
CSM Ocean Model, H is determined as the depth where the buoyancy difference 

g[/?(H)- /)(surface)]p(H)' 1 = 3 10' 4 ms' 2 (34/) 

XVII. Model constants 

As discussed in the Introduction, we use the new two— point closure turbulence model 
discussed in paper I to determine the critical time scales. The results are: 
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(35a) 


T pv r ' = 5' rr ptf= 5(1 + ^'). V' ,= ff t 

<r t =0.72, 7 i= i 0 = 0 . 6 , 

10a =1 +|f% a =6a , la =2-7 a , a =1 F=0.64 
5 5 i 5 ’ 2 2 5 3 5’ 


The values of the A's defined in Eq.(14d) are then as follows: 

Model A: (A ,A ,A ,A ,A ,A ,A ,A )= 

1 2 3 4 5 6 7 8 

(0.107, 0.0032, 0.0864, 0.12, 11.9, 0.4, 0, 0.48) 


(35b) 


The corresponding Ri cr is: 


Ri cr = 0.846 (35c) 

which is more than four times larger than the MY result. On the other hand, to lend 
further support to the turbulence model that led us to (35b), we have devised an 
alternative procedure (Appendix D) which is based on entirely different considerations and 
which gives very similar results: 

Model B: 


(A ,A ,A ,A ,A ,A ,A ,A ) — 

1 2 3 4 5 6 7 8 y 

(0.127, 0.00336, 0.0906, 0.101, 11.2, 0.4, 0, 0.318) 


with Ri : 
cr 


Ri cr =1.03 

In the first model, we have also used data from Shih and Shabbir (1992). 


(35d) 

(35e) 


XVIII. Tests of the model 

We begin with the case of pure shear, Ri-*0, which has been widely studied. In this 
case, we have (to first order) from (21a, c) 

DS m = V D = d o’ S m = 75 ( 36a ) 

Thus, the first of (20a) becomes 

V = c , t ’ c „= 011 (36b) 

which is the well known formula employed in shear flows studies (Rodi, 1984). 
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In Fig.l we exhibit the dimensionless structure functions S m ^ vs. Ri for the full local 
model, Eqs.(21a,c) and (26d). The left hand side of the figure is for unstable stratification 
(Ri<0) while the right hand panel is for stable stratification, Ri>0. As one can expect, the 
relative position of the curves switches as one moves from negative to positive Ri. In the 
unstable case, the heat diffusivity is larger than the momentum diffusivity since one 
expects that the temperature gradient affects more the heat diffusivity than the momentum 
diffusivity. The reverse is true when stratification is stable. The lower value of vs. S m is 
in accord with the laboratory data for the turbulent Prandtl number which we discuss 
below. In Fig. 2 we exhibit the dynamical time scale r in units of the shear E vs. Ri, that is 
the function ^’(Ri), Eq.(26c), which is solution of Eq.(26d). As expected, the stronger the 
stratification, the longer is the eddies life time, or equivalently, the smaller the turbulent 
kinetic energy, since In Fig. 3 we report the laboratory data for ct t vs. Ri by Webster 

(1964) and in Fig. 4 we exhibit the predicted turbulent Prandtl number a vs. Ri, Eqs.(20b) 
and (21). In Fig. 5 we show the predicted flux Richardson number R^ vs. Ri defined in 
Eq.(30c). In Fig. 6, we present a set of laboratory data for R^ vs. Ri from the work of 
Maderich et al. (1995). Since the definitions of the Richardson numbers in Figs. 5 and 6 are 
not identical, the two figures cannot be superimposed; however, what is important is the 
general behavior which is very similar. In Fig. 7 we exhibit the predicted value of the 
efficiency parameter 7 defined in Eq.(30a,c). 

In Fig. 8 we exhibit the rate of dissipation of turbulent kinetic energy € vs. depth z with 
and without the wave breaking phenomenon. As one can observe, the dotted line 
representing the law of the wall, Eq.(32a), provides a good representation of the data only 
when the wave breaking phenomenon is not important, which occurs at considerable 
depths. Near the surface, the data indicate a much larger dissipation rate (crosses and 
triangles, Terray et al., 1996) which the turbulence model can reproduce quite well (full 
line) using the dynamic equation for the turbulent kinetic energy (25) with the boundary 
condition at z=0 given by relation (32b). The values of H , u* and c are taken from tables 
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1,2 of Terray et al. ( 1996) and from table 1 of Drennan et al.( 1996). 


XIX. Ocean GCM 

To test the new vertical diffusivities, we used the global ocean general circulation 
model NCAR CSM Ocean Model produced by the University Corporation for Research, 
National Center for Atmospheric Research, Climate and Global Dynamics Division. They 
developed the model from the MOM 1.1 GFDL code (see the NCAR CSM Ocean Model 
Tech. Note, The NCAR CSM Ocean Model, by the NCAR Oceanography section). We 
employed the stand-alone 3°x3° configuration of the model detailed in their technical note 
with the default parameter values. It has 3.6° spacing in longitude and a variable spacing in 
latitude increasing from 1.8° at the equator to 3.4° at 17° N, S and then decreasing back to 
1.8° for 60° N, S and poleward. There are 25 levels of increasing thickness in the vertical, 
with the surface level 6 meters thick. The option for the GM mesoscale eddy 
parameterization was enabled. Bulk forcing with a seasonal cycle plus a 1/2 year timescale 
restoring condition on the salinity is used, except under sea— ice where there is strong 
restoring. This configuration corresponds to the case B— K described in Large et al. (1997). 
It should be noted, however, that for determination of the length scale in our turbulence 
model we used the program's definition for mixed layer depth (a buoyancy difference from 
the surface of 3 10' 4 ms' 2 ), which is different from that graphed as a diagnostic in Fig. 5 of 
Large et al.(1997). We initialized our runs with annually averaged Levitus data and ran for 
126 momentum years. As in Large et al. (1997) a 3504sec timestep for momentum is used, 
while for the first 96 momentum years the tracers are accelerated by a factor increasing 
from 10 at the surface to 100 for the deep ocean. We then set all timesteps equal for the 
remaining 30 years as they did. 

First, we ran the NCAR program as is, with the option for the KPP mixing enabled, 
producing the KPP data presented in the figures below. Then, in place of the KPP module, 
we inserted a module which uses our new model for the diffusivities for momentum and 
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heat with the salt diffusivity set equal to that of heat. To save computing time, we 
constructed tables of the dimensionless functions S m ^ and of the dimensionless variable y 
obtained from solving Eq.26a and using Eq.(25c), 

y = )x 2 (|) 2 = x 2 e A 2 £ 2 K‘ (37a) 


vs. Ri. Then, for each point in space and time these were interpolated to the local Ri. The 


diffusivities K m ^(model) were rewritten as 


Vt,/^ - * B ,r 5 s 


m,h 


(37c) 


XX. Below the Mixed Layer 

Below the ocean mixed layer, we employed the same background diffusivities as in the 
NCAR model, specifically: 

K =16.7 cmV\ K,=K =0.5 cmV 1 (38) 

m ’ h s 


XXI. Ocean GCM 

The results using the new vertical diffusivities K m and are presented in 

Figs. 9-20 where we exhibit the model results (squares) Levitus data (-) and the results 
using the KPP model (diamonds). We present both global T and S, Figs. 9-10, as well as 
specific basins, artic, atlantic, pacific, indian and southern oceans. 

XXII. Conclusions 

The goal of this paper was not only that of devising a model for the vertical 
diffusivities that performs better than previous models. The overall goal had a much wider 
breadth, that of building a model that satisfies several conditions: 

1) the model is not tailored specifically to ocean turbulence like previous model (e.g.,KPP), 

2) before being used in an O-GCM, it must be shown to reproduce well documented 
laboratory, atmospheric and LES data on stratified turbulence, e.g., Figs.3-4, 5-6, 
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3) it must incorporate the latest advances in turbulent closure modeling so as to naturally 
give rise to a critical Richardson number (Fig. 2) that corresponds not only to the empirical 
data from mixed layer studies (e.g., Martin, 1985) but which also reproduces laboratory 
data (e.g., Monin and Yaglom, 1971), 

4) in the non-local case, where the flux of turbulent kinetic energy is included in the 
equation for the kinetic energy, the model must be able to reproduce the observed increase 
of dissipation caused by breaking waves at the ocean surface, Fig. 8, 

5) the model not only reproduces the K^=7tN" 2 representation of the heat diffusivity that 
has been widely used in the past but it predicts a value for 7 that varies with Ri, as indeed 
expected, Fig. 7, 

6) after these tests have been passed, the model is tested in an O— GCM without changing 
any of the ingredients that have been used to reproduce 1)— 5) above. The model is expected 
to perform at least as well as previous ad hoc models like KPP. Even in the simplest case, 
corresponding to a local model, the results shown in Figs. 9— 20 satisfy this criterion, 

7) the strength of the model is its ability to encompass other cases within the same 
methodology. In the next paper, the model will be extended to include the salinity field 
thus allowing K g to be different than K^, as dictated by laboratory and ocean data. 
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Appendix A 

Pressure Correlation terms 


For the tensor 11-, we employ the most general expression (Canuto 1994): 


n.. = 2r- , b. ; + (1-/? )B- ■ - -KS-- — a E-. - o Z-. 
ij pv ij v V ij 5 ij i JJ 2 ij 


where we have defined the following traceless tensors and vectors: 

b.. = T..-I K«.., 

IJ ij 3 ij’ 


h i=¥ 


(A.l) 

(A. 2) 
(A. 3) 
(A.4) 

(A. 5) 

Sy = WU.J+UJ , ,). v if i(U iJ -U ji ,) (A. 6 ) 

In all the above equations, A-a^g. (g-=0,0,g), a T =-p-'dp/d T. Analogously, we have 
(Canuto, 1994) 


B ij“ A i h j + A j h i"^ij A k h k 
E ij ES ik b kj + S jk b ik ~ 3 ^ij S k/ b kf 

Z ij « b ik V b kj V ik 


n. = r;],0+7A/- a a (S-. + *V..)0 
i p» i 'i i 4 r ij 3 ij J 


(A. 7) 


Appendix B 

Third-order moments, TOM 

The equations for the third-order moments are taken from Canuto (1994). In the 
presence of buoyancy, shear and rotation, they are: 



= - (ujUjU^Uj^+perm) - (ujuJ^u^+perm)* (l-c^X^ftTjt^ + perm) 
-^ijq\ + perm ) (B-l) 

<St+Y'>W = yM - (¥k w j,k + Yk lW i,k> - 
<¥k i k w j + Yk i k TO i + w k i k ¥j> + 
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k A. A, ^u, + (1-c )(A.Fu-+A.^u.) 
3 a ij k k v n /v i J J v 


(B.2) 


( nt +r 3 ' 1+2r * 1) "I F = -" I v ; r i i ' r; Vij + 

+ (^ 11 )^-¥ j fe j F (B- 3 ) 

(Sr + = 3 ^/“j - m ik r + *-f> < B4 ) 

J «*/ 

where r =r/ 2c* and c*=7, c =4, c =1/5. In the absence of shear, and in the stationary 

3 10 11 

case, Eqs. (B.1)-(B.4) were solved analytically. The resulting third-order moments 
exhibited some unexpected symmetry properties in that they they are all given by the sum 
of the gradients of all the second— order moments, see Eq.(5d). These new expressions for 
the TOM were shown to reproduce quite well the LES data for a buoyancy driven PBL 
(Canuto et al., 1994). 


Appendix C 

Wave breaking case. Boundary condition for K 

Craig and Banner (1994) first suggested that at z=0 one takes the flux of turbulent 
kinetic energy F(K), see Eq.(7c), equal to 

F(K,z=0) = au*, a~10 2 (C.l) 

where u* is the wind stress at the surface. To translate (C.l) into a relation for the kinetic 
energy at z=0, we proceed in a different manner than Craig and Banner. Consider the 
stationary limit of Eq.(25) near the surface. In the absence of a heat flux, we obtain 

K m S2 + !< K m f > = f < C - 2 > 

Since N 2 ~0, Eq.(20a) for K m becomes K m = c KV 1 , c =2S m =constant and (C.2) becomes 

K m S2 + K nl< K m f > = C , K2 ( C 3 > 

Since near the surfacer K S 2 e - uwS = u 2 and since K 2 >>u 2 , Eq.(C.3) becomes 

m 13 
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a2 4 


d 2 K _ k2 i K d _ d 

Multiplying both sides of (C.4) by dK/d£, we obtain 

d ,dK, 2 _ 2 dK 3 

which can be integrated to yield 

K mff = AK /2 ’ 

Using the second of (C.l), we obtain 

ml = AK 3 / 2 

Since (Terray et al, 1996, Eq.5) 

ft = c/u* 

and using K=^q 2 , we finally obtain that at z=0 the velocity q is given by: 

q = (12c- 1 ) 1 / 6 u */ 3 c 1 / 3 
which is relation (32b) of the text. 


(C.4) 


(C.5) 

(C.6) 

(C.7) 

(C.8) 

(C.9) 


Appendix D 

Alternative Determination of the constants (14d) 

To determine the constants appearing in Eqs.(14e-f), we adopt from Shih and Shabbir 
(1992) and Canuto (1993) the expressions 

A =(l-a)(2cV, A =(l-a)(2cV 


A = 0 (2c Y\ 

4 5 4 


A = 1 - jQf , 

7 4 3 ’ 


A = 1 - |a 
6 4 3 

A =(1-7)7* 

8 1 T 


where 


(D.l) 


c* = 2 + 6.22 F 2 (l-F) 3 / 4 -F’ 

cn = 60 , Ot = i(2 — 7 a ) 

1 5 2 3 V 5 

a 5 = l0 ( 1+ | F? )> a 3 = t 

F = 0.64, 0 = 0.48 

5 

The values thus obtained are 


(D.2) 
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(A 2 , A , A , A , A ) = (0.00336, 0.0906, 0.101, 0.4, 0) (D.3) 

To determine A , we employ a neutral surface layer in which the mean profiles are 
1>5 

logarithmic, the mean wind is along the x-direction and z approaches zero, 

-uw = uj, 7-j^=(q/u*) 2 , iq 2 =K (D.4) 

Under these conditions, the algebraic equations for u 2 , v 2 and w 2 , Eqs.(18a— c), can be 
simplified and yield 

A = |(3A 2 - A 2 ) + 4B- 4 / 3 (D.5) 

where 

B 1 = (q/u*) 3 = 16.6 (D.6) 

in accord with the MY model. Using the values of A and A given in (D.3), we obtain 

A^ = 0.127 (D.7) 

The value of A is obtained in a similar fashion. Assuming that production equals 

5 

dissipation in the above neutral surface layer, we have 

Sh = 2B* /V: 1 , a. = 1 (D.8) 

1 to to 

where a is the turbulent Prandtl number in the neutral case which we take to be unitv. 
to 

Applying the neutral condition and z-*0 to the equations for u$, w$, Eqs.(19a,b), and using 

^=-K,S h dl (D.9) 

after some algebra we get an algebraic equation for A , 

A 2 - 7rB 4 / 3 (l + A -3 A )a. A - j(A -A )B 4 / 3 (A +A +2<r t ) = 0 (D.10) 

The solution is 

A = 11.2 (D.ll) 

Similarly, from the algebraic equation for IP, Eq.(8a), we have 

Y =g Q b j 2 / 3 Q, Qe^/(^ = 3.1 (D.12) 

where the subscript s indicates a surface quantity, and Q=3.1 in accordance with Mellor 
and Yamada (1982). Substituting (D.12) into (D.l) and using 7^1/3 gives 

A 8 = (i - 7j)^ = 0 318 (D. 13) 
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Figure caption 

Fig.l. Model predictions for the dimensionless structure functions S m and defined in 
Eq.(20a) vs. the degree of stability Ri. 

Fig. 2. The time scale r, specifically, the variable ip, Eq.(26c), vs. the degree of stability Ri, 
see Eq. (26c). Since Eq.(27) shows that ip is the inverse of the turbulent kinetic energy K, 
this figure shows that at Ri cr ~l, K vanishes for the stratification is to strong for turbulent 
mixing to survive. Several 1-point closure models of the Mellor— Yamada type predict a far 
smaller Ri cr “0.25 which considerably reduces the extent of turbulent mixing 
Fig. 3. Laboratory data for the turbulent Prandtl number vs.Ri (Webster 1964); a T is 
defined in Eq.(20b) and is the ratio of S m /S^ presented in Fig.l. 

Fig. 4. The turbulent Prandtl number a T = S m /S^ vs. Ri predicted by the present model. 
Several previous models predict a constant a ■ . 

Fig. 5. The predicted flux Richardson number R^ vs. Ri, see Eq.(30c). 

Fig. 6. The flux Richardson number Rr vs. Ri from a variety of laboratory data. The 
symbols refer to different authors cited in Maderich et al. (1995) and refer to either grid 
generated turbulence and/or freely decaying turbulence in a stably stratified medium. Since 
the Ri used in Figs. 5-6 are not identical, a superposition of the two figures in not possible. 
What is important is that the general behavior predicted by the model mirrors that of the 
data. 

Fig. 7. Efficiency parameter 7(Ri) defined in Eq.(30a,c) vs. Ri. Even though a great deal of 
caution must be exercised in using the data (see discussion after Eq.30c), an observed value 
of 7~l/4 (e.g., Oakey, 1982; Toole et al., 1994) is found to correspond to an Ri~0.25 which 
is well within the expected values. 

Fig. 8. The rate of dissipation of turbulent kinetic energy e(z) vs. depth z in meters. The 
predicted value of e(z) in the absence of the wave— breaking phenomenon (dashed line) is a 
poor representation of the data (crosses and diamonds) in the first 10 meters. A much more 
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realistic prediction (solid line) is achieved if wave-breaking is included in the turbulence 
model. This requires a non-local model, the details of which are discussed in Appendix C. 
Fig. 9. When the new expressions for K m ^ and K g =K^ are used in the O-GCM, we obtain 
the results exhibited in Figs. 9-20. To ease the comparison, we also ran the same O-GCM 
with the KPP model. The model results are compared with Levitus et al. (1994) data. 
Global temperature vs. depth. 

Fig. 10. Same as Fig. 9 for the global salinity 
Fig.l 1 . Same as Fig.9 for the Artie ocean 
Fig. 12. Same as Fig. 10 for the Artie ocean 
Fig. 13. Same as Fig.9 for the Atlantic ocean 
Fig. 14. Same as Fig. 10 for the Atlantic ocean 
Fig. 15. Same as Fig.9 for the Pacific ocean 
Fig. 16. Same as Fig. 10 for Pacific ocean 
Fig. 17. Same as Fig.9 for the Indian ocean 
Fig. 18. Same as Fig. 10 for the Indian ocean 
Fig. 19. Same as Fig.9 for the Southern ocean 
Fig. 20. Same as Fig. 10 for the Southern ocean 
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XIV. Passive Scalar 

In Appendix C we develop a model to compute the diffusivity K c for a passive scalar, 
e.g., tritium. The mean concentration C satisfies the equation 


DC _ d — 

DT ~ ” W u i c 


(30a) 


Here, u^c is the turbulence-induced concentration flux. In the stationary, local case, 


wc = - K J? 


where the passive scalar diffusivity K c is given by: 

K C = » J VP + (H,) V «^ 


(30b) 

(30c) 


with 


A = 1 + (IT, Vc^ (3M) 

Using Eqs.(18c) and (19a) for the turbulent variables w 2 and w D, we can express (30c) in 
the same functional form as (20a, b), namely 

K c = 2S c 7 ’ S c = 3X7 EC < 1 -^ A ) w 

with A>0 given by: 

A = A „ S m + A , S h Ri < 3 °8) 

A = 3A - A , A = 4A + 3(1— 7 )r a r A (30h) 

0 321 41 CV ' 

A = 1 + (l-\)r pc T ce r 2 m (30i) 

We recall that ^ is given by Eqs.(28c-e). As for the time scales, we suggest to identify 

T pc = V r 0c =r « and em P‘°y < 4f -e)- 

The case of a passive scalar 
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Following the formalism developed in Canuto (1994), the dynamic equation for the 


turbulent velocity u^ is given by: 


Du 

nr— Vu"' 


I. - jstVj _ u i u j> + V + 

1 J J J 

while that of the fluctuating temperature 6 is: 


d 2 u- 


D6 _ rp d , 0 — m , d 2 0 

m~~ u i T ,i _ 3x1V ~ V> + *t ^ 


(C.l) 


(C.2) 


If we further consider a passive scalar field C, the corresponding turbulent component c 
satisfies the equation: 


Dc _ n d , — \ . d 2 c 

nr - - "i c ,i - at j (u i c - u i c) + 


(C-3) 


where x c is the kinematic diffusivity of the c-field. To construct the required correlation 
wc, we multiply (C.l) by c, (C.3) by u- , average the results and sum the results. We 


obtain: 


where 


Dt cu i + D f = - + cu j u i,j) + - "i - £ c 

d 


D 


f = ^ cu i u j)’ n i s<d 3* i > ' f c Emi i,jj 


(C.4) 


(C.5) 


where we have used <..> instead of a bar for notational convenience and where e l 
represents the dissipation of cm due to molecular forces. For IlV we employ an expression 
analogous to (A. 7) but retain only the first two terms, 

n i = r P i ¥ + Vi^ (c - 6) 

Clearly, we need the correlation W, which we derive using Eqs.(C.2)— (C.3). The result is: 

(C.7) 


nt^ + D f - - H. + y $£.) + x T ^jj + A c ^jj 


The last two terms can be rearranged to read 

*T^,jj + = «X c + *T>TOjj - (X C + *T> + * 

(C.8) 

The first term may be incorporated with the diffusion term, the last term can be neglected 
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compared with others while the second term, taking into account that x T >>X g , can be 
parameterized as 




(C.9) 

so that finally 

+ D f = - <¥ H. + u j c fjj.) “ T c$ ^ 

(C.10) 

In the stationary, local (no diffusion) case, the model gives: 




(C. 11) 


r P J ¥ = - (¥j c ,j + + A i< 1 “ 7 ,) w 

(C.12) 

Thus, 

A u¥ = - ( K ctj ^ 

(C.13) 

where 

V V W u i,i + (1 -V A i r co T ,j> 

(C.14) 


< K c )ij E v [ Vj + (1 -V T c0 A i¥ l 

(C.15) 

In the case of vertical diffusivities, we further have 



— „ ac 

wc = - K cW 

(C. 16) 

where the passive scalar diffusivity K c is given by: 



K c = + ( 1 “9' 1 )t- c ^« x wJl 

(C.17) 

with 

A e 1 + (l- 7i ) Tpc r c N 2 

(C.18) 


We recall that w 2 and w# are provided by the turbulence model. 
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